home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / RLaB / test < prev    next >
Text File  |  1994-04-26  |  43KB  |  1,711 lines
  1. tic();
  2. //
  3. //  RLaB test program
  4. //
  5.  
  6. //
  7. // Test Parameters
  8. //
  9.  
  10. X = 3;
  11.  
  12. epsilon = function() 
  13. {
  14.   eps = 1.0;
  15.   while((eps + 1.0) != 1.0) 
  16.   {
  17.     eps = eps/2.0;
  18.   }
  19.   return eps;
  20. };
  21.  
  22. eps = epsilon();
  23.  
  24. eye = function( m , n ) 
  25. {
  26.   local(i, N, new);
  27.  
  28.   if (!exist (n))
  29.   {
  30.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  31.     new = zeros (m[1], m[2]);
  32.     N = min ([m[1], m[2]]);
  33.   else
  34.     if (class (m) == "string" || class (n) == "string") {
  35.       error ("eye(), string arguments not allowed");
  36.     }
  37.     if (max (size (m)) == 1 && max (size (n)) == 1)
  38.     {
  39.       new = zeros (m[1], n[1]);
  40.       N = min ([m[1], n[1]]);
  41.     else
  42.       error ("matrix arguments to eye() must be 1x1");
  43.     }
  44.   }
  45.   for(i in 1:N)
  46.   {
  47.     new[i;i] = 1.0;
  48.   }
  49.   return new;
  50. };
  51.  
  52. symm = function( A )
  53. {
  54.   return (A + A')./2;
  55. };
  56.  
  57. //
  58. //-------------------- Test Relational Expressions -------------------
  59. //
  60. printf("\tstart scalar tests...\n");
  61. printf("\tstart relational tests...\n");
  62.  
  63. //    SCALAR CONSTANTS (REAL)
  64. if( !(1<2) ) { error(); }
  65. if( !(1<=2) ) { error(); }
  66. if( 1>2 ) { error(); }
  67. if( 1>=2 ) { error(); }
  68. if( 1==2 ) { error(); }
  69. if( !(1!=2) ) { error(); }
  70. if( !1 ) { error(); }
  71. if( !!!1) { error(); }
  72.  
  73. if( !([1]<[2]) ) { error(); }
  74. if( !([1]<=[2]) ) { error(); }
  75. if( [1]>[2] ) { error(); }
  76. if( [1]>=[2] ) { error(); }
  77. if( [1]==[2] ) { error(); }
  78. if( !([1]!=[2]) ) { error(); }
  79. if( ![1] ) { error(); }
  80. if( !!![1]) { error(); }
  81.  
  82. //    SCALAR CONSTANTS (COMPLEX)
  83. if( !(1+2i<2+3i) ) { error(); }
  84. if( !(1+2i<=2+3i) ) { error(); }
  85. if( 1+2i>2+3i ) { error(); }
  86. if( 1+2i>=2+3i ) { error(); }
  87. if( 1+2i==2+3i ) { error(); }
  88. if( !(1+2i!=2+3i) ) { error(); }
  89. if( !1+2i ) { error(); }
  90. if( !!!1+2i) { error(); }
  91.  
  92. if( !([1+2i]<[2+3i]) ) { error(); }
  93. if( !([1+2i]<=[2+3i]) ) { error(); }
  94. if( [1+2i]>[2+3i] ) { error(); }
  95. if( [1+2i]>=[2+3i] ) { error(); }
  96. if( [1+2i]==[2+3i] ) { error(); }
  97. if( !([1+2i]!=[2+3i]) ) { error(); }
  98. if( ![1+2i] ) { error(); }
  99. if( !!![1+2i]) { error(); }
  100.  
  101. //    SCALAR ENTITIES (REAL)
  102. a=1;b=2;
  103. if( !(a<b) ) { error(); }
  104. if( !(a<=b) ) { error(); }
  105. if( a>b ) { error(); }
  106. if( a>=b ) { error(); }
  107. if( a==b ) { error(); }
  108. if( !(a!=b) ) { error(); }
  109. if( !a ) { error(); }
  110. if( !!!a) { error(); }
  111.  
  112. if( !([a]<[b]) ) { error(); }
  113. if( !([a]<=[b]) ) { error(); }
  114. if( [a]>[b] ) { error(); }
  115. if( [a]>=[b] ) { error(); }
  116. if( [a]==[b] ) { error(); }
  117. if( !([a]!=[b]) ) { error(); }
  118. if( ![a] ) { error(); }
  119. if( !!![a]) { error(); }
  120.  
  121. //    SCALAR ENTITIES (COMPLEX)
  122. a=1+2i;b=2+3i;
  123. if( !(a<b) ) { error(); }
  124. if( !(a<=b) ) { error(); }
  125. if( a>b ) { error(); }
  126. if( a>=b ) { error(); }
  127. if( a==b ) { error(); }
  128. if( !(a!=b) ) { error(); }
  129. if( !a ) { error(); }
  130. if( !!!a) { error(); }
  131.  
  132. if( !([a]<[b]) ) { error(); }
  133. if( !([a]<=[b]) ) { error(); }
  134. if( [a]>[b] ) { error(); }
  135. if( [a]>=[b] ) { error(); }
  136. if( [a]==[b] ) { error(); }
  137. if( !([a]!=[b]) ) { error(); }
  138. if( ![a] ) { error(); }
  139. if( !!![a]) { error(); }
  140.  
  141. if (! all (all (-2 <  rand(4,4)))) { error (); }
  142. if (! all (all (-2 <= rand(4,4)))) { error (); }
  143. if (! all (all ( 2 >  rand(4,4)))) { error (); }
  144. if (! all (all ( 2 >= rand(4,4)))) { error (); }
  145.  
  146. if (! all (all (rand(4,4) >  -2))) { error (); }
  147. if (! all (all (rand(4,4) >= -2))) { error (); }
  148. if (! all (all (rand(4,4) <   2))) { error (); }
  149. if (! all (all (rand(4,4) <=  2))) { error (); }
  150.  
  151. if (! all (all (rand (4,4) >  -rand (4,4)))) { error (); }
  152. if (! all (all (rand (4,4) >= -rand (4,4)))) { error (); }
  153. if (! all (all (-rand (4,4) <  rand (4,4)))) { error (); }
  154. if (! all (all (-rand (4,4) <= rand (4,4)))) { error (); }
  155.  
  156. //------------------------- Test REAL SCALARS ------------------------
  157. //
  158. //    CONSTANTS
  159. //      Addition
  160. if(1+2 != 3) { error(); }
  161. //      Subtraction
  162. if(1-2 != -1) { error(); }
  163. //      Multiply
  164. if(1*2 != 2) { error(); }
  165. //      Divide
  166. if(1/2 != 0.5) { error(); }
  167. //      Power
  168. if(2^2 != 4) { error(); }
  169. if(4^0 != 1) { error(); }
  170. //      Unary Minus
  171. if(2-3 != -1) { error(); }
  172. //
  173. //    ENTITIES
  174. //
  175. a = 1; b = 2; c = 3; d = 0.5;
  176. //      Addition
  177. if(a+b != c) { error(); }
  178. //      Subtraction
  179. if(a-b != -a) { error(); }
  180. //      Multiply
  181. if(a*b != b) { error(); }
  182. //      Divide
  183. if(a/b != d) { error(); }
  184. //      Power
  185. if(b^b != b*b) { error(); }
  186. if(b^0 != 1) { error (); }
  187. //      Unary Minus
  188. if(b-c != -a) { error(); }
  189. //
  190. //  ENTITIES & CONSTANTS
  191. //
  192. if(a+2 != c) { error(); }
  193. //      Subtraction
  194. if(1-b != -a) { error(); }
  195. //      Multiply
  196. if(a*2 != 2) { error(); }
  197. //      Divide
  198. if(1/b != d) { error(); }
  199. //      Power
  200. if(2^b != b*b) { error(); }
  201. //      Unary Minus
  202. if(b-3 != -a) { error(); }
  203. //
  204. //------------------------Test COMPLEX SCALARS -------------------------
  205. //
  206. //    CONSTANTS
  207. if(sqrt(-1) != 1i) { error(); }
  208. //      Addition
  209. if((1+2i)+(2+3i) != (3+5i)) { error(); }
  210. //      Subtraction
  211. if((1+2i)-(3+4i) != (-2-2i)) { error(); }
  212. //      Multiply
  213. if((1+2i)*(3+4i) != (-5+10i)) { error(); }
  214. //      Divide
  215. if((1+2i)/(3-4i) != (-.2+.4i)) { error(); }
  216. //      Power
  217. //      Precision problems prevent us from testing these. Have to
  218. //      be checked by hand.
  219. //  (1+2i)^2 = -3 + 4i
  220. //  (1+2i)^.5 = 1.272 + 7.862e-1i
  221. //  if((1+2i)^2 != (-3+4i)) { error(); }
  222. //  if((1+2i)^10 != (237+3116i)) { error(); }
  223. //      Unary Minus
  224. if(-(1+2i) != -1-2i) { error(); }
  225. //
  226. //    ENTITIES
  227. //
  228. a = 1+2i; b = 3+4i; c = 4+6i;
  229. if(a+b != c) { error(); }
  230. //      Subtraction
  231. if(a-b != -2-2i) { error(); }
  232. //      Multiply
  233. if(a*b != -5+10i) { error(); }
  234. //      Divide
  235. //if(a/(3-4i) != -.2+.4i) { error(); }
  236. //      Power
  237. //  if(b^b != b*b) { error(); }
  238. //      Unary Minus
  239. if(-a != -1-2i) { error(); }
  240. //
  241. //    ENTITIES & CONSTANTS
  242. //
  243. if(a+(3+4i) != c) { error(); }
  244. //      Subtraction
  245. if((1+2i)-b != -2-2i) { error(); }
  246. //      Multiply
  247. if(a*(3+4i) != -5+10i) { error(); }
  248. //      Divide
  249. //if(a/(3-4i) != -.2+.4i) { error(); }
  250. //      Power
  251. //if(b^b != b*b) { error(); }
  252. //      Unary Minus
  253. if(-(1+2i) != -1-2i) { error(); }
  254. //
  255. // String - Numerical Equalities
  256. //
  257. if ((1 == "1")) { error(); }
  258. if (([1] == "1")) { error(); }
  259. if (("1" == 1)) { error(); }
  260. if (("1" == [1])) { error(); }
  261.  
  262. if (rand(3,3) == "str") { error(); }
  263. if ("str" == rand(3,3)) { error(); }
  264. if (!any (any ((["1", "2"; "3", "4"] == "1") == [1,0;0,0]))) { error(); }
  265.  
  266. //
  267. //----------------------- Test REAL MATRICES ---------------------------
  268. //
  269. printf("\tstart matrix tests...\n");
  270. printf("\t\treal-matrices\n");
  271. //  Read in test matrices
  272. //
  273. read("test.input");
  274. //
  275. //  Matrix construction
  276. //
  277. if(all([3] != matrix(3))) { 
  278.   error(); 
  279. }
  280. if(any([1;2;3] != [1,2,3]')) {
  281.   error();
  282. }
  283. if(any (any (m0 != zeros(2,2)))) { error(); }
  284. if(any (any (m1 != 1+zeros(2,2)))) { error(); }
  285. if(any (any (m2 != [1,2;3,4]))) { error(); }
  286. if(any (any (m3 != [1+2i,2+3i;3+4i,5+6i]))) { error(); }
  287. if(any (any ([1,2;3+0i,4+0i] != m2))) { error(); }
  288. if(any (any ([m2] != m2))) { error(); }
  289. //
  290. //  Sub-Matrix promotion
  291. //
  292. if(m2[2;2] != 4) { error(); }
  293. if(any(m2[2;] != [3,4])) { error(); }
  294. if(any(m2[;2] != [2,4]')) { error(); }
  295. i=2;j=1;
  296. if(m2[i;j] != 3) { error(); }
  297. i=1;j=2;
  298. if(m2[i;j] != 2) { error(); }
  299. m = [1,2,3;4,5,6;7,8,9];
  300. if(any(m[1;1,2] != [1,2])) {
  301.   error();
  302. }
  303. if(any(m[1,2;1] != [1;4])) {
  304.   error();
  305. }
  306. if(any (any (m[1,2;1,2] != [1,2;4,5]))) {
  307.   error();
  308. }
  309. //
  310. if(m3[2;2] != (5+6i)) { error(); }
  311. if(any(m3[2;] != [3+4i,5+6i])) { error(); }
  312. if(any(m3[;2] != conj([2+3i,5+6i]'))) { error(); }
  313. //
  314. //  Automatic creation, extension
  315. //
  316. if(any (any ((mm[3;3]=10) != [0,0,0;0,0,0;0,0,10]))) { error(); }
  317. a=[1,2,3;4,5,6];
  318. if(any (any ((a[3;1]=10) != [1,2,3;4,5,6;10,0,0]))) { error(); }
  319. a=[1,2;3,4];
  320. if(any (any ((a[3,4;3,4]=[5,6;7,8]) != [1,2,0,0;3,4,0,0;0,0,5,6;0,0,7,8]))) {
  321.   error();
  322. }
  323. mmm[2;] = a[1;];
  324. //
  325. //  Matrix binary operations
  326. //
  327. a = m2; b = [5,6;7,8];
  328. if(any (any (a+a != [2,4;6,8]))) { error(); }
  329. if(any (any (a-a != zeros(2,2)))) { error(); }
  330. if(any (any (2+a != [3,4;5,6]))) { error(); }
  331. if(any (any (2-a != [1,0;-1,-2]))) { error(); }
  332. if(any (any (a-2 != [-1,0;1,2]))) { error(); }
  333. if(any (any (2*a != [2,4;6,8]))) { error(); }
  334. if(any (any ((a./2 != [0.5,1;1.5,2])))) { error(); }
  335. if(any (any (a*a != [7,10;15,22]))) { error(); }
  336. if(any (any (a*a*a != [37,54;81,118]))) { error(); }
  337. if(any (any (a .* a != [1,4;9,16]))) { error(); }
  338. if(any (any (a./a != [1,1;1,1]))) { error(); }
  339. if(any (any (a' != [1,3;2,4]))) { error(); }
  340.  
  341. if(any(any(rand(3,3)^0 != eye(3,3)))) { error(); }
  342. if(any(any(rand(3,3).^0 != ones(3,3)))) { error(); }
  343. if(any(any(rand(1,3).^0 != ones(1,3)))) { error(); }
  344. if(any(any(rand(3,1).^0 != ones(3,1)))) { error(); }
  345. if(any(any(1.^zeros(3,1) != ones(3,1)))) { error(); }
  346. if(any(any(1.^zeros(1,3) != ones(1,3)))) { error(); }
  347.  
  348. if(any ([1;2;3]+[4;5;6] != [5;7;9])) {
  349.   error();
  350. }
  351. if(any ([1;2;3]-[4;5;6] != [-3;-3;-3])) {
  352.   error();
  353. }
  354. if(any ([2;2;2] ./ [1;1;1] != [2;2;2])) {
  355.   error();
  356. }
  357. if(any ([1;2;3] .* [4;5;6] != [4;10;18])) {
  358.   error();
  359. }
  360.  
  361. //
  362. // Test row-wise matrix addition
  363. //
  364.  
  365. a = [1,2,3]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  366. if (!all (all (a .+ b == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  367. if (!all (all (b .+ a == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  368.  
  369. a = [1,2,3] + [1,2,3]*1i; 
  370. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  371. c = [2,4,6;5,7,9;8,10,12;11,13,15] + [2,4,6;5,7,9;8,10,12;11,13,15]*1i;
  372. if (!all (all (a .+ b == c))) { error (); }
  373. if (!all (all (b .+ a == c))) { error (); }
  374.  
  375. printf("\t\tpassed matrix row-wise add test...\n");
  376.  
  377. //
  378. // Test row-wise matrix subtraction
  379. //
  380.  
  381. a = [1,1,1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  382. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  383. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  384.  
  385. a = [1,1,1] + [1,1,1]*1i;
  386. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  387. c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
  388. if (!all (all (a .- b == -c))) { error (); }
  389. if (!all (all (b .- a ==  c))) { error (); }
  390.  
  391. printf("\t\tpassed matrix row-wise subtraction test...\n");
  392.  
  393. //
  394. // Test col-wise matrix addition
  395. //
  396.  
  397. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  398. if (!all (all (a .+ b == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  399. if (!all (all (b .+ a == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  400.  
  401. a = [1;1;1;1] + [1;1;1;1]*1i;
  402. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  403. c = [2,3,4;5,6,7;8,9,10;11,12,13] + [2,3,4;5,6,7;8,9,10;11,12,13]*1i;
  404. if (!all (all (a .+ b == c))) { error (); }
  405. if (!all (all (b .+ a == c))) { error (); }
  406.  
  407. printf("\t\tpassed matrix col-wise add test...\n");
  408.  
  409. //
  410. // Test col-wise matrix subtraction
  411. //
  412.  
  413. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  414. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  415. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  416.  
  417. a = [1;1;1;1] + [1;1;1;1]*1i;
  418. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  419. c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
  420. if (!all (all (a .- b == -c))) { error (); }
  421. if (!all (all (b .- a ==  c))) { error (); }
  422.  
  423. printf("\t\tpassed matrix col-wise subtraction test...\n");
  424.  
  425. a = [1,2,3];
  426. b = [1,2,3;4,5,6;7,8,9];
  427. c = [1,4,9;4,10,18;7,16,27];
  428.  
  429. if (!all (all (a .* b == c))) { error (); }
  430. if (!all (all (b .* a == c))) { error (); }
  431.  
  432. za = a + rand (size (a))*1j;
  433. zb = b + rand (size (b))*1j;
  434.  
  435. if (!all (all (za .* zb == [za;za;za] .* zb))) { error (); }
  436. if (!all (all (zb .* za == zb .* [za;za;za]))) { error (); }
  437. if (!all (all (a .* zb == [a;a;a] .* zb))) { error (); }
  438. if (!all (all (zb .* a == zb .* [a;a;a]))) { error (); }
  439. if (!all (all (za .* b == [za;za;za] .* b))) { error (); }
  440. if (!all (all (b .* za == b .* [za;za;za]))) { error (); }
  441.  
  442. printf("\t\tpassed matrix row-wise multiplication test...\n");
  443.  
  444. a = [1,2,3];
  445. b = [1,2,3;4,6,6;7,8,9];
  446. c = [1,1,1;4,3,2;7,4,3];
  447.  
  448. if (!all (all (b ./ a == c))) { error (); }
  449. if (!all (all ([a;a;a] ./ b == a ./ b))) { error (); }
  450. if (!all (all (b ./ [a;a;a] == b ./ a))) { error (); }
  451.  
  452. za = a + rand (size (a))*1j;
  453. zb = b + rand (size (b))*1j;
  454.  
  455. if (!all (all ([za;za;za] ./ zb == za ./ zb))) { error (); }
  456. if (!all (all (zb ./ [za;za;za] == zb ./ za))) { error (); }
  457. if (!all (all ([a;a;a] ./ zb == a ./ zb))) { error (); }
  458. if (!all (all (zb ./ [a;a;a] == zb ./ a))) { error (); }
  459. if (!all (all ([za;za;za] ./ b == za ./ b))) { error (); }
  460. if (!all (all (b ./ [za;za;za] == b ./ za))) { error (); }
  461.  
  462. printf("\t\tpassed matrix row-wise division test...\n");
  463.  
  464. a = [1;2;3];
  465. b = [1,2,3;4,5,6;7,8,9];
  466.  
  467. if (!all (all (a .* b == [a,a,a] .* b))) { error (); }
  468. if (!all (all (b .* a == b .* [a,a,a]))) { error (); }
  469.  
  470. za = a + rand (size (a))*1j;
  471. zb = b + rand (size (b))*1j;
  472.  
  473. if (!all (all (za .* zb == [za,za,za] .* zb))) { error (); }
  474. if (!all (all (zb .* za == zb .* [za,za,za]))) { error (); }
  475. if (!all (all (za .* b == [za,za,za] .* b))) { error (); }
  476. if (!all (all (b .* za == b .* [za,za,za]))) { error (); }
  477. if (!all (all (a .* zb == [a,a,a] .* zb))) { error (); }
  478. if (!all (all (zb .* a == zb .* [a,a,a]))) { error (); }
  479.  
  480. printf("\t\tpassed matrix column-wise multiplication test...\n");
  481.  
  482. a = [1;2;3];
  483. b = [1,2,3;4,6,6;7,8,9];
  484.  
  485. if (!all (all ([a,a,a] ./ b == a ./ b))) { error (); }
  486. if (!all (all (b ./ [a,a,a] == b ./ a))) { error (); }
  487.  
  488. za = a + rand (size(a))*1j;
  489. zb = b + rand (size(b))*1j;
  490.  
  491. if (!all (all ([za,za,za] ./ zb == za ./ zb))) { error (); }
  492. if (!all (all (zb ./ [za,za,za] == zb ./ za))) { error (); }
  493. if (!all (all ([za,za,za] ./ b == za ./ b))) { error (); }
  494. if (!all (all (b ./ [za,za,za] == b ./ za))) { error (); }
  495. if (!all (all ([a,a,a] ./ zb == a ./ zb))) { error (); }
  496. if (!all (all (zb ./ [a,a,a] == zb ./ a))) { error (); }
  497.  
  498. printf("\t\tpassed matrix column-wise division test...\n");
  499.  
  500.  
  501. //
  502. //--------------------- Test COMPLEX MATRICES --------------------------
  503. //
  504. //  Automatic creation, extension
  505. //
  506. printf("\t\tcomplex-matrices\n");
  507. if(any (any ((mm[3;3]=10+10i) != [0,0,0;0,0,0;0,0,10+10i]))) { error(); }
  508. a=[1,2,3;4,5,6];
  509. if(any (any ((a[3;1]=10+10i) != [1,2,3;4,5,6;10+10i,0,0]))) { error(); }
  510. //
  511. a = m3;
  512. if(any (any (a+a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  513. if(any (any (a-a != zeros(2,2)))) { error(); }
  514. if(any (any (2+a != [3+2i,4+3i;5+4i,7+6i]))) { error(); }
  515. if(any (any (2-a != [1-2i,0-3i;-1-4i,-3-6i]))) { error(); }
  516. if(any (any (a-2 != [-1+2i,0+3i;1+4i,3+6i]))) { error(); }
  517. if(any (any (2*a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  518. if(any (any (a./2 != [.5+1i,1+1.5i;1.5+2i,2.5+3i]))) { error(); }
  519. if(any (any (a*a != [-9+21i,-12+34i;-14+48i,-17+77i]))) { error(); }
  520. if(any (any (a*a*a != [-223+57i,-345+113i;-469+183i,-719+337i]))) { error(); }
  521. if(any (any (a .* a != [-3+4i,-5+12i;-7+24i,-11+60i]))) { error(); }
  522. //
  523. // The following test may not work on some machines
  524. //
  525. if(any (any (a./a != [1,1;1,1]))) { 
  526.   printf("\t\t***complex division inaccuracy, check manually***\n");
  527. }
  528.  
  529. if(any (any (a' != conj([1+2i,3+4i;2+3i,5+6i])))) { error(); }
  530. //  
  531. //--------------------- Test NULL MATRICES -------------------------
  532. //
  533. printf("\t\tnull-matrices\n");
  534. // Create a NULL matrix
  535. mnull = matrix();
  536. // Test it with SCALARS
  537. if( any([1,mnull] != 1)) {
  538.   error();
  539. }
  540. if( any([mnull,1] != 1)) {
  541.   error();
  542. }
  543. // Test with MATRIX construction
  544. m = [1,2;3,4;5,6];
  545. if( any([mnull;1] != [1])) {
  546.   error();
  547. }
  548. if( any([1;mnull] != [1])) {
  549.   error();
  550. }
  551. if( any([mnull;1,2,3] != [1,2,3])) {
  552.   error();
  553. }
  554. if( any([1,2,3;mnull] != [1,2,3])) {
  555.   error();
  556. }
  557. if(any (any ([mnull,m] != m))) {
  558.   error();
  559. }
  560. if(any (any ([m,mnull] != m))) {
  561.   error();
  562. }
  563. if(any (any ([mnull;m] != m))) {
  564.   error();
  565. }
  566. if(any (any ([m;mnull] != m))) {
  567.   error();
  568.  
  569. mnull = matrix();
  570. mnull[1] = [1];
  571. }
  572.  
  573. //--------------------- Test Matrix Multiply --------------------------
  574.  
  575. i = sqrt(-1);
  576. a = [1,2,3;4,5,6;7,8,9];
  577. b = [4,5,6;7,8,9;10,11,12];
  578. c = [ 48,  54,  60 ;
  579.      111, 126, 141 ;
  580.      174, 198, 222 ];
  581.  
  582. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  583.  
  584. az = a + b*i;
  585. bz = b + a*i;
  586.  
  587. cz = [-18+141*i , -27+162*i , -36+183*i ;
  588.         9+240*i ,   0+279*i ,  -9+318*i ;
  589.        36+339*i ,  27+396*i ,  18+453*i ];
  590.  
  591. czz = [ 48+30*i ,  54+36*i  ,  60+42*i ;
  592.        111+66*i , 126+81*i  , 141+96*i ;
  593.        174+102*i, 198+126*i , 222+150*i ];
  594.  
  595. czzz = [ 48+111*i ,  54+126*i ,  60+141*i ;
  596.         111+174*i , 126+198*i , 141+222*i ;
  597.         174+237*i , 198+270*i , 222+303*i ];
  598.  
  599. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  600. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  601. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  602.  
  603. a = [a,a];
  604. b = [b;b];
  605. c = [  96 , 108 , 120 ;
  606.       222 , 252 , 282 ;
  607.       348 , 396 , 444 ];
  608.  
  609. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  610.  
  611. az = [az,az];
  612. bz = [bz;bz];
  613.  
  614. cz = [  -36+282*i ,  -54+324*i ,  -72+366*i ;
  615.          18+480*i ,    0+558*i ,  -18+636*i ;
  616.          72+678*i ,   54+792*i ,   36+906*i ];
  617.  
  618. czz = [  96+60*i  , 108+72*i  , 120+84*i  ;
  619.         222+132*i , 252+162*i , 282+192*i ;
  620.         348+204*i , 396+252*i , 444+300*i ];
  621.  
  622. czzz = [  96+222*i , 108+252*i , 120+282*i ;
  623.          222+348*i , 252+396*i , 282+444*i ;
  624.          348+474*i , 396+540*i , 444+606*i ];
  625.  
  626. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  627. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  628. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  629.  
  630. printf("\t\tpassed matrix multiply test...\n");
  631.  
  632. //--------------------- Test STRING MATRICES --------------------------
  633. //
  634. printf("\t\tstring-matrices\n");
  635. sm = ["s1","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"];
  636. if(sm[1] != "s1") { error(); }
  637. if( sm[1;3] != "sm3" ) { error(); }
  638. if(any(sm[2,3;3] != ["xxx3";"yyy3"]) ) { error(); }
  639. if(any (any ((sm[1;1]="xx")!=["xx","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"])))
  640. {
  641.   error();
  642. }
  643. if( ((strm[1] = "strm")[1]) != "strm" ) { error(); }
  644.  
  645. //  Test string-matrix add functionality
  646.  
  647. sm = sm[1,2;1,2];
  648. if (any (any (("1_"+sm+"_2") != ["1_xx_2","1_sm2_2";"1_x1_2","1_x2_2"]))) {error();}
  649.  
  650. //
  651. //---------------------------- List Tests --------------------------
  652. //
  653. //  List creation
  654. listest = << << 11; 12 >>; << 21; 22>> >>;
  655. if( listest.[1].[2] != 12 ) { error(); }
  656. if(any(<<a=10;b=1:4;c=[1,2,3;4,5,6;7,8,9]>>.b != [1,2,3,4])) { error(); }
  657. mlist.[0] = m;
  658. if(any(any(mlist.[0] != m))) { error(); }
  659. printf("\tlist tests\n");
  660.  
  661. //
  662. // Reset random number generator seed...
  663. //
  664.  
  665. rand("default");
  666.  
  667. //
  668. //-------------------------- Test abs () ------------------------------
  669. //
  670.  
  671. A = rand(5,5);
  672. T = ( A == abs (A));
  673. if (!all (all (A))) { error ("abs() incorrect"); }
  674. printf("\tabs()");
  675.  
  676. //
  677. //-------------------------- Test max () ------------------------------
  678. //
  679.  
  680. A = [1,10,100;2,20,200;1,2,3];
  681. if (!all (max (A) == [2,20,200])) { error( "max() incorrect"); }
  682. if (max (max(A)) != 200) { error ("max() incorrect"); }
  683. printf("\tmax()");
  684.  
  685. //
  686. //-------------------------- Test min () ------------------------------
  687. //
  688.  
  689. if (!all (min (A) == [1,2,3])) { error( "min() incorrect"); }
  690. if (min (min(A)) != 1) { error ("min() incorrect"); }
  691. printf("\tmin()");
  692.  
  693. //
  694. //-------------------------- Test maxi () -----------------------------
  695. //
  696.  
  697. if (!all (maxi (A) == [2,2,2])) { error( "maxi() incorrect"); }
  698. if (maxi (maxi(A)) != 1) { error ("maxi() incorrect"); }
  699. printf("\tmaxi()");
  700.  
  701. //
  702. //-------------------------- Test mini () -----------------------------
  703. //
  704.  
  705. if (!all (mini (A) == [1,3,3])) { error( "mini() incorrect"); }
  706. if (mini (mini(A)) != 1) { error ("mini() incorrect"); }
  707. printf("\tmini()");
  708.  
  709. //
  710. //-------------------------- Test floor () ----------------------------
  711. //
  712.  
  713. if (floor (1.9999) != 1) { error ("floor() output incorrect"); }
  714. if (!all (all (floor ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  715.   error ("floor output incorrect");
  716. }
  717. printf("\tfloor()");
  718.  
  719. //
  720. //-------------------------- Test ceil () 0----------------------------
  721. //
  722.  
  723. if (ceil (1.9999) != 2) { error ("ceil() output incorrect"); }
  724. if (!all (all (ceil ([1.99,1.99;2.99,2.99]) == [2,2;3,3]))) {
  725.   error ("ceil output incorrect");
  726. }
  727. printf("\tceil()\n");
  728.  
  729. //
  730. //-------------------------- Test round () ----------------------------
  731. //
  732.  
  733. if (round (1.8) != 2) { error ("round() output incorrect"); }
  734. if (round (1.4) != 1) { error ("round() output incorrect"); }
  735. if (!all (all (round ([1.99,1.99;2.4,2.4]) == [2,2;2,2]))) {
  736.   error ("round output incorrect");
  737. }
  738. printf("\tround()");
  739.  
  740. //
  741. //-------------------------- Test sum () ------------------------------
  742. //
  743.  
  744. S = [1:4; 4:7; 8:11];
  745. if (sum (S[1;]) != 10) { error ("sum() incorrect"); }
  746. if (sum (S[3;]) != 38) { error ("sum() incorrect"); }
  747. if (!all (all (sum (S) == [13,16,19,22]))) { error ("sum() incorrect"); }
  748. printf("\tsum()");
  749.  
  750. //
  751. //-------------------------- Test int () ------------------------------
  752. //
  753.  
  754. if (int (1.9999) != 1) { error ("int() output incorrect"); }
  755. if (!all (all (int ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  756.   error ("int() output incorrect");
  757. }
  758. printf("\tint()");
  759.  
  760. //
  761. //-------------------------- Test mod () ------------------------------
  762. //
  763.  
  764. if (mod (1,1) != 0) { error ("mod() output incorrect"); }
  765. if (mod (4,2) != 0) { error ("mod() output incorrect"); }
  766. if (mod (3,2) != 1) { error ("mod() output incorrect"); }
  767. if (mod (5,3) != 2) { error ("mod() output incorrect"); }
  768. printf("\tmod()\n");
  769.  
  770. //
  771. //-------------------------- Test rand () -----------------------------
  772. //
  773. rand ("normal", 5, 1);
  774. xrand = rand(4000,1);
  775.  
  776. mean = function(x)
  777. {
  778.   local(m);
  779.  
  780.   m = size (x)[1];
  781.   if( m == 1 ) 
  782.   { 
  783.     m = size (x)[2];
  784.   }
  785.  
  786.   return sum( x ) / m;
  787. };
  788.  
  789. std = function(x)
  790. {
  791.   local(i, m, s);
  792.  
  793.   if(class(x) != "num") { error("std() requires NUMERICAL input"); }
  794.  
  795.   m = x.nr;
  796.   if( m == 1 ) 
  797.   { 
  798.     return sqrt( sum( (x - mean(x)) .^ 2 ) / (x.nc - 1) );
  799.   else
  800.     for( i in 1:x.nc) {
  801.       s[i] = sqrt( sum( (x[;i] - mean(x[;i])) .^ 2 ) / (x.nr - 1) );
  802.     }
  803.     return s;
  804.   }
  805. };
  806.  
  807. if (!(mean (xrand) > 4.9 && mean (xrand) < 5.1)) 
  808.   { error ("error in random"); }
  809. if (!(std (xrand) > 0.9 && std (xrand) < 1.1))
  810.   { error ("error in random"); }
  811. printf("\tpassed rand test...\n");
  812.  
  813. rand("default");
  814.  
  815. //
  816. //-------------------------- Test norm () -----------------------------
  817. //
  818.  
  819. tn = [1,2,3,4;2,1,2,3;3,2,1,2;4,3,2,1  ];
  820. if (norm(tn,"m") != 4 ) { error ("incorrect norm computation"); }
  821. if (norm(tn,"1") != 10) { error ("incorrect norm computation"); }
  822. if (norm(tn,"i") != 10) { error ("incorrect norm computation"); }
  823. printf("\tpassed norm test...\n");
  824.  
  825. //
  826. //-------------------------- Test qr () ------------------------------
  827. //
  828.  
  829. a = rand (4, 4);
  830. qa = qr (a);
  831. if (max (max (abs (qa.q*qa.r - a)))/(X*norm (a)*a.nr) > eps)
  832.   { error ("possible qr() problems"); }
  833.  
  834. z = rand (4,4) + rand(4,4)*1i;
  835. qz = qr (z);
  836. if (max (max (abs (qz.q*qz.r -  z)))/(X*norm (z)*z.nr) > eps)
  837.   { error ("possible qr() problems"); }
  838.  
  839. printf("\tpassed qr test...\n");
  840.  
  841. //
  842. //-------------------------- Test schur () ----------------------------
  843. //
  844.  
  845. a = rand (4, 4);
  846. sa = schur (a);
  847. if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
  848.   { error ("possible schur() problems"); }
  849.  
  850. z = rand (4,4) + rand(4,4)*1i;
  851. sz = schur (z);
  852. if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
  853.   { error ("possible schur() problems"); }
  854.  
  855. printf("\tpassed schur test...\n");
  856.  
  857. //
  858. //-------------------------- Test schord () ----------------------------
  859. //
  860.  
  861. s2a = schord (sa, 2, 4);
  862. if (max (max (abs (s2a.z*s2a.t*s2a.z' - a)))/(X*norm (a)*a.nr) > eps)
  863.   { error ("possible schord() problems"); }
  864.  
  865. s2z = schord (sz, 3, 1);
  866. if (max (max (abs (s2z.z*s2z.t*s2z.z' - z)))/(X*norm (z)*z.nr) > eps)
  867.   { error ("possible schord() problems"); }
  868.  
  869. printf("\tpassed schord test...\n");
  870.  
  871.  
  872. //
  873. //-------------------------- Test chol () -----------------------------
  874. //
  875.  
  876. c = rand(4,4);
  877. c = symm (c + eye (size (c))*norm(c)*c.nr);
  878. u = chol (c);
  879. if (max (max (abs (u'*u - c)))/(X*norm (c)*c.nr) > eps)
  880.   { error ("possible chol() problems"); }
  881.  
  882. cz = rand(4,4);
  883. cz = symm (cz + eye (size (cz))*X*norm(cz)*cz.nr);
  884. uz = chol (cz);
  885. if (max (max (abs (uz'*uz - cz)))/(X*norm (cz)*cz.nr) > eps)
  886.   { error ("possible chol() problems"); }
  887.  
  888. printf("\tpassed chol test...\n");
  889.  
  890. //
  891. //-------------------------- Test solve () -----------------------------
  892. //
  893.  
  894. a = rand(4,4);
  895. while (1/rcond (a) > 100) {  a = rand(4,4); }
  896. b = ones(4,1);
  897. x = solve (a,b);
  898. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  899.   { 
  900.     printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  901.     printf ("\tA*X - B:\n");
  902.     abs (a*x - b)
  903.     error ("possible solve() problems\n");
  904.   }
  905.  
  906. az = rand(4,4) + rand(4,4)*1j;
  907. while (1/rcond (az) > 100) { az = rand(4,4) + rand(4,4)*1j; }
  908. bz = rand(4,1) + rand(4,1)*1j;
  909. xz = solve (az,bz);
  910. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  911.   { 
  912.     printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  913.     printf ("\tA*X - B:\n");
  914.     abs (az*xz - bz)
  915.     error ("possible solve() problems\n");
  916.   }
  917.  
  918. printf("\tpassed solve test...\n");
  919.  
  920.  
  921. //
  922. //-------------------------- Test lu() / factor() -----------------------------
  923. //
  924.  
  925. tril = function(x, k) 
  926. {
  927.   local(i, j, nr, nc, y);
  928.  
  929.   if (!exist (k)) { k = 0; }
  930.   nr = x.nr; nc = x.nc;
  931.   if(k > 0) 
  932.   { 
  933.     if (k > (nc - 1)) { error ("tril: invalid value for k"); }
  934.   else
  935.     if (abs (k) > (nr - 1)) { error ("tril: invalid value for k"); }
  936.   }
  937.  
  938.   y = zeros(nr, nc);
  939.  
  940.   for(i in max( [1,1-k] ):nr) {
  941.     j = 1:min( [nc, i+k] );
  942.     y[i;j] = x[i;j];
  943.   }
  944.  
  945.   return y;
  946. };
  947.  
  948. triu = function(x, k) 
  949. {
  950.   local(i, j, nr, nc, y);
  951.  
  952.   if (!exist (k)) { k = 0; }
  953.   nr = x.nr; nc = x.nc;
  954.  
  955.   if(k > 0) 
  956.   { 
  957.     if (k > (nc - 1)) { error ("triu: invalid value for k"); }
  958.   else
  959.     if (abs (k) > (nr - 1)) { error ("triu: invalid value for k"); }
  960.   }
  961.  
  962.   y = zeros(nr, nc);
  963.  
  964.   for(j in max( [1,1+k] ):nc) {
  965.     i = 1:min( [nr, j-k] );
  966.     y[i;j] = x[i;j];
  967.   }
  968.  
  969.   return y;
  970. };
  971.  
  972. static (swap);
  973.  
  974. lu = function ( A )
  975. {
  976.   local (i, l, u, pvt, x)
  977.  
  978.   if (A.nr != A.nc) { error ("lu() requires square A"); }
  979.  
  980.   x = factor (A);    // Do the factorization
  981.  
  982.   //
  983.   // Now create l, u, and pvt from lu and pvt.
  984.   //
  985.  
  986.   l = tril (x.lu, -1) + eye (size (x.lu));
  987.   u = triu (x.lu);
  988.   pvt = eye (size (x.lu));
  989.  
  990.   //
  991.   // Now re-arange the columns of pvt
  992.   //
  993.  
  994.   for (i in 1:max (size (x.lu)))
  995.   {
  996.     pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
  997.   }
  998.   return << l = l; u = u; pvt = pvt >>;
  999. };
  1000.  
  1001. //
  1002. //  In vector V, swap elements I, J
  1003. //
  1004.  
  1005. swap = function ( V, I, J )
  1006. {
  1007.   local (v, tmp);
  1008.   v = V;
  1009.   tmp = v[I];
  1010.   v[I] = v[J];
  1011.   v[J] = tmp;
  1012.   return v;
  1013. };
  1014.  
  1015.  
  1016. a = rand(4,4);
  1017. while (1/rcond (a) > 100) {  a = rand(4,4); }
  1018. lua = lu (a);
  1019. if (max (max (abs(a - lua.pvt*lua.l*lua.u)))/(X*norm (a)*a.nr) > eps)
  1020.   { 
  1021.     printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  1022.     printf ("\tA - p*l*u:\n");
  1023.     abs (a - lua.pvt*lua.l*lua.u)
  1024.     error ("possible lu()/factor() problems\n");
  1025.   }
  1026.  
  1027. az = rand(4,4) + rand(4,4)*1j;
  1028. while (1/rcond (az) > 100) { az = rand(4,4) + rand(4,4)*1j; }
  1029. luz = lu (az);
  1030. if (max (max (abs (az - luz.pvt*luz.l*luz.u)))/(X*norm (az)*az.nr) > eps)
  1031.   { 
  1032.     printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  1033.     printf ("\tA - p*l*u:\n");
  1034.     abs (az - luz.ovt*luz.l*luz.u)
  1035.     error ("possible lu()/factor()() problems\n");
  1036.   }
  1037.  
  1038. printf("\tpassed lu/factor test...\n");
  1039.  
  1040.  
  1041. //
  1042. //-------------------------- Test svd ()   -----------------------------
  1043. //
  1044.  
  1045. a = rand(4,4);
  1046. s = svd (a);
  1047. if (max (max (abs (s.u*diag(s.sigma)*s.vt - a)))/(X*norm (a)*a.nr) > eps)
  1048.   { error ("possible svd() problems"); }
  1049.  
  1050. z = rand(4,4) + rand(4,4)*1j;
  1051. sz = svd (z);
  1052. if (max (max (abs (sz.u*diag(sz.sigma)*sz.vt - z)))/(X*norm (z)*z.nr) > eps)
  1053.   { error ("possible svd() problems"); }
  1054.  
  1055. printf("\tpassed svd test...\n");
  1056.  
  1057.  
  1058. //
  1059. //-------------------------- Test hess ()  -----------------------------
  1060. //
  1061.  
  1062. a = rand(4,4);
  1063. h = hess (a);
  1064. if (max (max (abs (h.p*h.h*h.p' - a)))/(X*norm (a)*a.nr) > eps)
  1065.   { error ("possible hess() problems"); }
  1066.  
  1067. z = rand(4,4) + rand(4,4)*1j;
  1068. hz = hess (z);
  1069. if (max (max (abs (hz.p*hz.h*hz.p' - z)))/(X*norm (z)*z.nr) > eps)
  1070.   { error ("possible hess() problems"); }
  1071.  
  1072. printf("\tpassed hess test...\n");
  1073.  
  1074. //
  1075. //-------------------------- Test lyap () ------------------------------
  1076. //
  1077.  
  1078. lyap = function ( A, B, C )
  1079. {
  1080.   local (X, sa, sb, tc)
  1081.  
  1082.   if (!exist (B)) 
  1083.   { 
  1084.     B = A';    // Solve the special form: A*X + X*A' = -C
  1085.   }
  1086.  
  1087.   if ((A.nr != A.nc) || (B.nr != B.nc) || (C.nr != A.nr) || (C.nc != B.nr)) {
  1088.     error ("Dimensions do not agree.");
  1089.   }
  1090.  
  1091.   //
  1092.   // Schur decomposition on A and B
  1093.   //
  1094.  
  1095.   sa = schur (A);
  1096.   sb = schur (B);
  1097.  
  1098.   //
  1099.   // transform C
  1100.   //
  1101.  
  1102.   tc = sa.z' * C * sb.z;
  1103.  
  1104.   X = sylv (sa.t, sb.t, tc);
  1105.  
  1106.   //
  1107.   // Undo the transformation
  1108.   //
  1109.  
  1110.   X = sa.z * X * sb.z';
  1111.  
  1112.   return X;
  1113. };
  1114.  
  1115. a = rand (5,5);
  1116. b = rand (5,5);
  1117. c = rand (5,5);
  1118. while (1/rcond (a) > 100) {  a = rand(5,5); }
  1119. while (1/rcond (b) > 100) {  b = rand(5,5); }
  1120. while (1/rcond (c) > 100) {  c = rand(5,5); }
  1121.  
  1122. x = lyap (a, b, c);
  1123. if (max (max (abs (a*x + x*b + c)))/(X*norm(c)*norm(a)*norm(b)) > eps)
  1124.   { error ("possible problems with lyap() or sylv()"); }
  1125.  
  1126. printf("\tpassed lyap test...\n");
  1127.  
  1128. //
  1129. //-------------------------- Test eig () ------------------------------
  1130. //
  1131.  
  1132. trace = function(m) 
  1133. {
  1134.   local(i, tr);
  1135.  
  1136.   if(m.class != "num") { 
  1137.     error("must provide NUMERICAL input to trace()");
  1138.   }
  1139.  
  1140.   tr = 0;
  1141.   for(i in 1:min( [m.nr, m.nc] )) {
  1142.     tr = tr + m[i;i];
  1143.   }
  1144.  
  1145.   return tr;
  1146. };
  1147.  
  1148. eye = function( m , n ) 
  1149. {
  1150.   local(i, N, new);
  1151.  
  1152.   if (!exist (n))
  1153.   {
  1154.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  1155.     new = zeros (m[1], m[2]);
  1156.     N = min ([m[1], m[2]]);
  1157.   else
  1158.     if (class (m) == "string" || class (n) == "string") {
  1159.       error ("eye(), string arguments not allowed");
  1160.     }
  1161.     if (max (size (m)) == 1 && max (size (n)) == 1)
  1162.     {
  1163.       new = zeros (m[1], n[1]);
  1164.       N = min ([m[1], n[1]]);
  1165.     else
  1166.       error ("matrix arguments to eye() must be 1x1");
  1167.     }
  1168.   }
  1169.   for(i in 1:N)
  1170.   {
  1171.     new[i;i] = 1.0;
  1172.   }
  1173.   return new;
  1174. };
  1175.  
  1176. //
  1177. // Standard eigenvalue problem
  1178. //
  1179.  
  1180. a = rand(10,10);
  1181. sa = symm (a);
  1182. ta = trace (a);
  1183. tsa = trace (sa);
  1184.  
  1185. z = rand(10,10) + rand(10,10)*1i;
  1186. sz = symm (z);
  1187. tz = trace (z);
  1188. tsz = trace (sz);
  1189.  
  1190. tol = 1.e-7;
  1191.  
  1192. if (!(ta < sum(eig(a).val) + tol && ta > sum(eig(a).val) - tol))
  1193. {
  1194.   error ("error in eig");
  1195. }
  1196. if (!(tsa < sum(eig(sa).val) + tol && tsa > sum(eig(sa).val) - tol))
  1197. {
  1198.   error ("error in eig");
  1199. }
  1200. if (!(tz < sum(eig(z).val) + tol && tz > sum(eig(z).val) - tol))
  1201. {
  1202.   error ("error in eig");
  1203. }
  1204. if (!(tsz < sum(eig(sz).val) + tol && tsz > sum(eig(sz).val) - tol))
  1205. {
  1206.   error ("error in eig");
  1207. }
  1208.  
  1209. //
  1210. // Generalized eigenvalue problem
  1211. //
  1212.  
  1213. b = rand(10,10);
  1214. sb = symm (b) + eye(size(b))*3;
  1215. tb = trace (b);
  1216. tsb = trace (sb);
  1217.  
  1218. zb = rand(10,10) + rand(10,10)*1i;
  1219. szb = symm (zb) + eye(size(zb))*3;
  1220. tzb = trace (zb);
  1221. tszb = trace (szb);
  1222.  
  1223. eig(a,b);    // not sure of a good way to check these yet
  1224. eigs(sa,sb);
  1225. eigs(sa,sb);
  1226.  
  1227. eig(z, zb);
  1228. eigs(sz, szb);
  1229. eigs(sz, szb);
  1230.  
  1231. printf("\tpassed eig test...\n");
  1232.  
  1233. //
  1234. //-------------------------- Test fft () -----------------------------
  1235. //
  1236.  
  1237. if (100 != fft(ones(100,1))[1]) { error ("error in fft()"); }
  1238. printf("\tpassed fft test...\n");
  1239.  
  1240. //
  1241. //-------------------------- Test printf () --------------------------
  1242. //
  1243.  
  1244. sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", 5,3,2, 8,7,3, "string", 3, 4, 1234e-2);
  1245. if (!(tmp == "  002  0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
  1246. printf("\tpassed sprintf test...\n");
  1247.  
  1248. //
  1249. //------------------------- Test strtod()  ----------------------------
  1250. //
  1251.  
  1252. if (123.456 != strtod ("123.456")) { error (); }
  1253. if (!all (all ([1,2;3,4] == strtod (["1","2";"3","4"]))))
  1254.   { error (); }
  1255. printf("\tpassed strtod test...\n");
  1256.  
  1257. //
  1258. //------------------------- Test getline()  ---------------------------
  1259. //
  1260. //
  1261.  
  1262. close( "test.getline" );
  1263.  
  1264. x = getline( "test.getline" );
  1265. if ( x.[1] !=  123.456 ) { error(); }
  1266. if ( x.[2] != -123.456 ) { error(); }
  1267. if ( x.[3] !=  123.456 ) { error(); }
  1268.  
  1269. x = getline( "test.getline" );
  1270. if ( x.[1] !=  .123 ) { error(); }
  1271. if ( x.[2] != -.123 ) { error(); }
  1272. if ( x.[3] !=  .123 ) { error(); }
  1273.  
  1274. x = getline( "test.getline" );
  1275. if ( x.[1] !=  123 ) { error(); }
  1276. if ( x.[2] != -123 ) { error(); }
  1277. if ( x.[3] !=  123 ) { error(); }
  1278.  
  1279. x = getline( "test.getline" );
  1280. if ( x.[1] !=  1e6 ) { error(); }
  1281. if ( x.[2] != -1e6 ) { error(); }
  1282. if ( x.[3] !=  1e6 ) { error(); }
  1283.  
  1284. x = getline( "test.getline" );
  1285. if ( x.[1] !=  1e5 ) { error(); }
  1286. if ( x.[2] != -1e5 ) { error(); }
  1287. if ( x.[3] !=  1e5 ) { error(); }
  1288.  
  1289. x = getline( "test.getline" );
  1290. if ( x.[1] !=  123.456e3 ) { error(); }
  1291. if ( x.[2] != -123.456e3 ) { error(); }
  1292. if ( x.[3] !=  123.456e3 ) { error(); }
  1293.  
  1294. x = getline( "test.getline" );
  1295. if ( x.[1] !=  123.456e3 ) { error(); }
  1296. if ( x.[2] != -123.456e3 ) { error(); }
  1297. if ( x.[3] !=  123.456e3 ) { error(); }
  1298.  
  1299. x = getline( "test.getline" );
  1300. if ( x.[1] !=  123.456e-3 ) { error(); }
  1301. if ( x.[2] != -123.456e-3 ) { error(); }
  1302. if ( x.[3] !=  123.456e-3 ) { error(); }
  1303.  
  1304. x = getline( "test.getline" );
  1305. if ( x.[1] !=  .123e3 ) { error(); }
  1306. if ( x.[2] != -.123e3 ) { error(); }
  1307. if ( x.[3] !=  .123e3 ) { error(); }
  1308.  
  1309. x = getline( "test.getline" );
  1310. if ( x.[1] !=  123e3 ) { error(); }
  1311. if ( x.[2] != -123e3 ) { error(); }
  1312. if ( x.[3] !=  123e3 ) { error(); }
  1313.  
  1314. x = getline( "test.getline" );
  1315. if ( x.[1] != "123abc" ) { error(); }
  1316. if ( x.[2] != "abc123" ) { error(); }
  1317. if ( x.[3] != "123.abc" ) { error(); }
  1318.  
  1319. x = getline( "test.getline" );
  1320. if ( x.[1] != "quoted string" ) { error(); }
  1321. if ( x.[2] != "q string with escapes \n \t \" " ) { error(); }
  1322.  
  1323. x = getline( "test.getline" );
  1324. if ( x.[1] != "quoted string" ) { error(); }
  1325. if ( x.[2] !=  1.23e3 ) { error(); }
  1326. if ( x.[3] !=  100 ) { error(); }
  1327. if ( x.[4] != "q string with escapes \n \t \" " ) { error(); }
  1328. if ( x.[5] !=  200 ) { error(); }
  1329.  
  1330. printf("\tpassed getline() test...\n");
  1331.  
  1332. //
  1333. //---------------------- Test readb()/writeb()  --------------------
  1334. //
  1335. a = rand (5,5);
  1336. z = rand(3,3) + rand(3,3)*1j;
  1337. strm = what()[1:5;1:5];
  1338. pi = 4*atan(1.0);
  1339. sc = 2*pi;
  1340. str = "this is a sample string\ttab";
  1341. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1342.  
  1343. writeb ("jnk.rb", a, z, strm, sc, str, l);
  1344.  
  1345. # Set aside matrices for later tests
  1346.  
  1347. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1348.  
  1349. clear (a, z, strm, sc, str, l);
  1350.  
  1351. close ("jnk.rb");
  1352.  
  1353. readb ("jnk.rb");
  1354.  
  1355. #
  1356. # Now do checks
  1357. #
  1358.  
  1359. if (!all (all (a == check.a))) { error (); }
  1360. if (!all (all (z == check.z))) { error (); }
  1361. if (!all (all (strm == check.strm))) { error (); }
  1362. if (sc != check.sc) { error (); }
  1363. if (str != check.str) { error (); }
  1364.  
  1365. if (length (l) != 5) { error (); }
  1366. if (!all (all (l.a == check.a))) { error (); }
  1367. if (!all (all (l.z == check.z))) { error (); }
  1368. if (!all (all (l.strm == check.strm))) { error (); }
  1369. if (l.sc != check.sc) { error (); }
  1370. if (l.str != check.str) { error (); }
  1371.  
  1372.  
  1373. printf("\tpassed binary I/O test...\n");
  1374.  
  1375. //
  1376. //---------------------- Test read()/write()  --------------------
  1377. //
  1378. a = rand (5,5);
  1379. z = rand(3,3) + rand(3,3)*1j;
  1380. strm = what()[1:5;1:5];
  1381. pi = 4*atan(1.0);
  1382. sc = 2*pi;
  1383. str = "this is a sample string\ttab";
  1384. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1385.  
  1386. write ("jnk.ra", a, z, strm, sc, str, l);
  1387.  
  1388. # Set aside matrices for later tests
  1389.  
  1390. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1391.  
  1392. clear (a, z, strm, sc, str, l);
  1393.  
  1394. close ("jnk.ra");
  1395.  
  1396. read ("jnk.ra");
  1397.  
  1398. #
  1399. # Now do checks
  1400. #
  1401.  
  1402. if (a.nr != 5 || a.nc != 5) { error (); }
  1403. if (z.nr != 3 || z.nc != 3) { error (); }
  1404. if (strm.nr != 5 || strm.nc != 5) { error (); }
  1405. if (str != check.str) { error (); }
  1406.  
  1407. if (length (l) != 5) { error (); }
  1408. if (l.a.nr != 5 || l.a.nc != 5) { error (); }
  1409. if (l.z.nr != 3 || l.z.nc != 3) { error (); }
  1410. if (l.strm.nr != 5 || l.strm.nc != 5) { error (); }
  1411. if (l.str != check.str) { error (); }
  1412.  
  1413. printf("\tpassed ascii I/O test...\n");
  1414.  
  1415. //
  1416. //------------------------- Fibonacci Test -------------------------
  1417. //
  1418. //  Calculate Fibonacci numbers
  1419. //
  1420. i=1; 
  1421. while ( i < 2 ) { 
  1422.   i=i+1;
  1423.   a=0; b=1;
  1424.   while ( b < 10000 ) {
  1425.     c = b;
  1426.     b = a+b;
  1427.     a = c;
  1428.   }
  1429. }
  1430. if ( b != 10946 ) {
  1431.   error("failed fibonacci test");
  1432. else
  1433.   printf("\tpassed fibonacci test...\n");
  1434. }
  1435.  
  1436. //
  1437. //------------------------- Factorial Test -------------------------
  1438. //
  1439. fac = function(a) {
  1440.   if(a <= 1) {
  1441.     return 1
  1442.   else
  1443.     return a*fac(a-1)
  1444.   }
  1445. };
  1446. if(fac(10) != 3628800) { error(); else printf("\tpassed factorial test...\n"); }
  1447. //
  1448. //--------------------------- ACK Test ----------------------------
  1449. //
  1450. ack = function(a, b) {
  1451.   if(a == 0) { return b + 1; }
  1452.   if(b == 0) { return $self(a-1, 1); }
  1453.   return $self(a-1, $self(a, b-1));
  1454. };
  1455.  
  1456. if(ack(2,2) != 7) {
  1457.   error("error in ack() test");
  1458. else
  1459.   printf("\tpassed ACK test...\n");
  1460. }
  1461. //
  1462. //------------------------- Prime Test -----------------------------
  1463. //
  1464. // An example that finds all primes less than limit
  1465. //
  1466. primes = function (limit) {
  1467.   local(prime, cnt, i, j, k);
  1468.  
  1469.   i = 1; cnt = 0;
  1470.   for(k in 2:limit) {
  1471.     j = 2;
  1472.     while(mod(k,j) != 0) {
  1473.       j++;
  1474.     }
  1475.     if(j == k) {            // Found prime
  1476.       cnt++;
  1477.       prime[i;1] = k;
  1478.       i++;
  1479.     }
  1480.   }
  1481.   return prime;
  1482. };
  1483.  
  1484. if(max(size(primes(100))) != 25) { 
  1485.   error("error in prime test");
  1486. else
  1487.   printf("\tpassed prime test...\n");
  1488. }
  1489.  
  1490. //
  1491. //--------------------------- Fib Min Test -----------------------------
  1492. //
  1493. //    fibmin() will minimize an arbitrary function 
  1494. //    in 1D using Fibonacci search
  1495.  
  1496. f065 = function(x)
  1497. {
  1498.     return (x - 0.65) * (x - 0.65);
  1499. };
  1500.  
  1501. fib = function(x)
  1502. {
  1503.     local (n, a, b);
  1504.  
  1505.     a = 1;
  1506.     b = 1;
  1507.     if (x >= 2)
  1508.     {
  1509.         n = x - 1;
  1510.         for (n in n:1:-1)
  1511.         {
  1512.             c = b;
  1513.             b = a + b;
  1514.             a = c;
  1515.             n = n - 1;
  1516.         }
  1517.     }
  1518.     return b;
  1519. };
  1520.  
  1521. //  Minimize a 1D function using Fibonacci search
  1522. //  f = function to minimize
  1523. //  xlo = lower bound
  1524. //  xhi = upper bound
  1525. //  n = number of iterations (the bigger the more accurate)
  1526. fibmin = function(f, xlo, xhi, n) {
  1527.     local(a, b, x, y, ex, ey, k, lo, hi);
  1528.  
  1529.     lo = xlo;
  1530.     hi = xhi;
  1531.     k = n;
  1532.     for (k in k:2:-1)
  1533.     {
  1534.         a = fib(k - 2) / fib(k);
  1535.                 b = fib(k - 1) / fib(k);
  1536.                 x = lo + (hi - lo) * a;
  1537.                 y = lo + (hi - lo) * b;
  1538.                 ex = f(x);
  1539.                 ey = f(y);
  1540.                 if (ex >= ey)
  1541.                 {
  1542.                         lo = x;
  1543.                 else
  1544.                         hi = y;
  1545.                 }
  1546. //  printf("%d: (%g %g) %g %g %g %g\n",  k, a, b, lo, hi, ex, ey);
  1547.         }
  1548.     return (lo + hi) / 2;
  1549. };
  1550.  
  1551. //
  1552. // Simple example using above function to mimize f065. Answer is 0.65
  1553. //
  1554. x = fibmin(f065, 0, 1, 30); // printf("f(%g)=%g\n", x, f065(x));
  1555. if (abs(x - 0.65) > 1e-6)
  1556. {
  1557.     printf("x = %f\n", x);
  1558.     error("failed fibmin test");
  1559. }
  1560.  
  1561. printf("\tpassed fibmin test...\n");
  1562.  
  1563. //
  1564. //--------------------- Nasty Function Test ------------------------
  1565. //
  1566. printf("\tStarting Nasty Function Test...");
  1567. printf("\tthis will take awhile\n");
  1568. check = function( a, b, c, d, e, f, g, h ) {
  1569.   if ( a+b+c+d == e+f+g+h && ...
  1570.        a^2+b^2+c^2+d^2 == e^2+f^2+g^2+h^2 && ...
  1571.        a^3+b^3+c^3+d^3 == e^3+f^3+g^3+h^3 ) {
  1572.     return 1;
  1573.   else
  1574.     return 0;
  1575.   }
  1576. };
  1577.  
  1578. for(a in 8:10) {
  1579.   for(b in 7:(a-1)) {
  1580.     for(c in 6:(b-1)) {
  1581.       for(d in 5:(c-1)) {
  1582.         for(e in 4:(d-1)) {
  1583.           for(f in 3:(e-1)) {
  1584.             for(g in 2:(f-1)) {
  1585.               for(h in 1:(g-1)) {                  
  1586.               if(check( a, b, c, d,  e, f, g, h ) || ...
  1587.                      check( a, e, c, d,  b, f, g, h ) || ...
  1588.                      check( a, f, c, d,  e, b, g, h ) || ...
  1589.                      check( a, g, c, d,  e, f, b, h ) || ...
  1590.                      check( a, h, c, d,  e, f, g, b ) || ...
  1591.                      check( a, b, e, d,  c, f, g, h ) || ...
  1592.                      check( a, b, f, d,  e, c, g, h ) || ...
  1593.                      check( a, b, g, d,  e, f, c, h ) || ...
  1594.                      check( a, b, h, d,  e, f, g, c ) || ...
  1595.                      check( a, b, c, e,  d, f, g, h ) || ...
  1596.                      check( a, b, c, f,  e, d, g, h ) || ...
  1597.                      check( a, b, c, g,  e, f, d, h ) || ...
  1598.                      check( a, b, c, h,  e, f, g, d ) || ...
  1599.                      check( a, e, f, d,  b, c, g, h ) || ...
  1600.                      check( a, e, g, d,  b, f, c, h ) || ...
  1601.                      check( a, e, h, d,  b, f, g, c ) || ...
  1602.                      check( a, f, g, d,  e, b, c, h ) || ...
  1603.                      check( a, f, h, d,  e, b, g, c ) || ...
  1604.                      check( a, g, h, d,  e, f, b, c ) || ...
  1605.                      check( a, b, e, f,  c, d, g, h ) || ...
  1606.                      check( a, b, e, g,  c, f, d, h ) || ...
  1607.                      check( a, b, e, h,  c, f, g, d ) || ...
  1608.                      check( a, b, f, g,  e, c, d, h ) || ...
  1609.                      check( a, b, f, h,  e, c, g, d ) || ...
  1610.                      check( a, b, g, h,  e, f, c, d ) || ...
  1611.                      check( a, e, f, g,  e, f, g, h ) || ...
  1612.                      check( a, e, f, h,  e, f, g, h ) || ...
  1613.                      check( a, e, g, h,  e, f, g, h ) || ...
  1614.                      check( a, f, g, h,  e, f, g, h ) ) { cnt++; }
  1615.               }
  1616.             }
  1617.           }
  1618.         }
  1619.       }
  1620.     }
  1621.   }
  1622. }
  1623.  
  1624. if(1) {  // figure out the value of cnt, and check!
  1625.   printf("\tpassed nasty function test...\n");
  1626. else
  1627.   error();
  1628. }
  1629.  
  1630. //
  1631. //------------------ Test More Advanced Functions --------------------
  1632. //
  1633.  
  1634. printf( "\tStarting the lqr/ode test..." );
  1635. printf( "\tthis will take awhile\n" );
  1636.  
  1637. lqr = function( a, b, q, r ) {
  1638.     local( k, s,...
  1639.            m, n, mb, nb, mq, nq,...
  1640.            e, v, d );
  1641.  
  1642.     m = size(a)[1]; n = size(a)[2];
  1643.     mb = size(b)[1]; nb = size(b)[2];
  1644.     mq = size(q)[1]; nq = size(q)[2];
  1645.     
  1646.     if ( m != mq || n != nq ) {
  1647.     fprintf( "stderr", "A and Q must be the same size.\n" );
  1648.     quit
  1649.     }
  1650.     
  1651.     mr = size(r)[1]; nr = size(r)[2];
  1652.     if ( mr != nr || nb != mr ) {
  1653.     fprintf( "stderr", "B and R must be consistent.\n" );
  1654.     quit
  1655.     }
  1656.     
  1657.     nn = zeros( m, nb );
  1658.     
  1659.     // Start eigenvector decomposition by finding eigenvectors of Hamiltonian:
  1660.     
  1661.     e = eig( [ a, solve(r',b')'*b'; q, -a' ] );
  1662.     v = e.vec; d = e.val;
  1663.     
  1664.     index = sort( real( d ) ).ind;
  1665.     d = real( d[ index ] );
  1666.     
  1667.     if ( !( d[n] < 0 && d[n+1] > 0 ) ) {
  1668.     fprintf( "stderr", "Can't order eigenvalues.\n" );
  1669.     quit
  1670.     }
  1671.     
  1672.     chi = v[ 1:n; index[1:n] ];
  1673.     lambda = v[ (n+1):(2*n); index[1:n] ];
  1674.     s = -real(solve(chi',lambda')');
  1675.     k = solve( r, nn'+b'*s );
  1676.     
  1677.     return << k=k; s=s >>;
  1678.  
  1679. };
  1680.  
  1681. // Now run a little test problem.
  1682.  
  1683. k = 1; m = 1; c = .1;
  1684. a = [0     ,1    ,0    , 0;
  1685.     -k/m, -c/m,  k/m,  c/m;
  1686.      0,     0,    0,    1;
  1687.      k/m,  c/m, -k/m, -c/m ];
  1688. b = [ 0; 1/m; 0; 0 ];
  1689. qxx = diag( [0, 0, 100, 0] );
  1690. ruu = [1];
  1691. K = lqr( a, b, qxx, ruu ).k;
  1692.  
  1693. dot = function( t, x ) 
  1694. {
  1695.     return (a-b*K)*x + b*K*([1,0,1,0]');
  1696. };
  1697.  
  1698. x = ode ( dot, [0,0,0,0], 0, 15, .02, 1e-5, 1e-5 );
  1699.  
  1700. m = maxi( x[;2] );
  1701.  
  1702. if ( (abs( x[m;2] - 1.195 ) > 0.001)  || ...
  1703.      any (abs( x[x.nr;2,4] - 1 ) > 0.001) ) {
  1704.      printf( "\tfailed***\n" );
  1705. else
  1706.      printf( "\tpassed the lqr/ode test...\n" );
  1707. }
  1708.  
  1709. printf("Elapsed time = %i seconds\n", toc() );
  1710. "FINISHED TESTS"
  1711.